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1.   Introduction 

The  finite  element  method  (FEM) ,  which  was  developed  in  engineering 
statics,  has  recently  been  introduced  into  various  atmospheric  prediction 
models  (Cullen,  197*1;  Hinsman,  1975;  Staniforth  and  Mitchell,  1977).  The 
FEM  is  a  special  case  of  the  Galerkin  procedure  in  which  the  dependent 
variables  are  approximated  by  a  finite  sum  of  spatially  varying  basis 
functions  with  time  dependent  coefficients.   The  FEM  basis  functions  are 
low  order  polynomials  which  are  zero  except  in  a  localized  region.   The 
Galerkin  procedure  produces  a  set  of  coupled  ordinary  differential 
equations  for  the  coefficients  which  are  solved  by  introducing  finite 
differences  in  time  (see  for  example  Pinder  and  Gray  (1977)). 

FEM  models  are  potentially  more  accurate  than  finite  difference 
models,  but  they  normally  require  more  computational  effort  per  degree  of 
freedom.   For  this  reason  it  is  especially  important  to  formulate  FEM 
models  efficiently.   Kelley  and  Williams  (1976)  found  considerable  small 
scale  noise  in  an  FEM  model  of  the  shallow  water  equations  in  a  channel 
which  had  all  variables  carried  at  the  same  nodal  points.  Winninghoff 
(1968),  Arakawa  and  Lamb  (1977)  and  Schoenstadt  (1980)  have  demonstrated 
the  advantages  of  spatial  staggering  of  dependent  variables  in  finite 
difference  models.  Also  Staniforth  and  Mitchell  (1977,  1978)  have 
obtained  excellent  results  with  a  vorticity-divergence  FEM  formulation. 
This  paper  will  compare  these  FEM  formulations  by  considering  the 
geostrophic  adjustment  process  with  the  linearized  shallow  water 
equations  in  one  dimension. 


2.   Basic  Equations 

The  linearized  shallow-water  equations  with  no  mean  flow  can  be 
written : 

|i+fu.O,  <2-2> 

t3t 

^  +  H^=0,  (2-3) 

where  u  and  v  are  the  perturbation  velocities  in  the  x  and  y  directions, 
respectively,  and  H  and  h  the  mean  and  perturbed  heights  of  the  free 
surface.   Also  g  represents  gravity  and  f  is  the  coriolis  parameter.   Note 
that  all  quantities  are  independent  of  y. 

The  vorticity  and  divergence  equations  are  obtained  by  differentiating 
(2.1)  and  (2.2)  with  respect  to  x  which  yields: 

|E.fc  +  gl!|=o,  (2.4) 

dx 

S+fD-O,  (2-5) 

3t 

2±+   HD  =  0  ,  (2-6) 

at 

where  D  =  3u/3x  is  the  divergence  and  X>   =  3v/3x  is  the  vorticity.   These 
relations  for  D  and  L,   are  particularly  simple  in  this  case  since  3u/3y  = 
3v/3y  =  0. 

Schoenstadt  (1977)  solved  the  continuous  equations  (2.1)-(2.3)  with 
the  of  the  spatial  Fourier  transform.   If  we  denote  Fourier  transforms  by 
a  tilde,  such  as 


u(k,t)  =   /     u(x,t)  e  ikx  dx  ,  (2.7) 


then  the  set  (2.1)-(2.3)  can  be  transformed  to  the  form: 

#  =  nfv  -  iygh  ,  (2-8) 

dt 

f=-nfS,  (2.9) 

4^=  -iyHu  ,  (2-10) 

dt 

where  H  =  1  and  y  =  k.   The  quantities  n  and  y  will  be  useful  later  when 
finite  difference  and  finite  element  solutions  are  needed.   The  initial 
conditions  are  written 


u     =  u(k,0)  =   /     u(x,0)  e  lkx  dx  ,         (2.11) 


with  similar  definitions  for  v  and  h  .   Schoenstadt  (1977)  solved  the  set 

o  o 

(2.7)-(2.9)    subject  to   initial   conditions  by  the  eigenvalue-eigenvector 
approach  which  gives : 

v                         ipgho 
u(k,t)    =  u     cos  Vt  +  nf  —  sin  vt   -  — sin  vt    ,  (2.12) 

'  O  V  v 

2  2   2 

v(k,t)    =  -  ^  u     sin  vt  +  {^-f-  +  Z-f-  cos  vt}   vo 

v        °  V  V 

+  iSmf.  {i   _   cos   Vt}   hQ    ,  (2.13) 

v 

h(k,t)    --M5      sin  vt   -  ^  {1   -   cos   vt}v 
'  v        o  v2  o 

+{H_f_  +  iL_|H   cos  vt}   hQ    ,  (2.14) 

V  v 


where:  2  2 .2  2  TT  o   -\r\ 

v=nf+ygH.  (2. IS; 


The  transformed  vorticity-divergence  set  (2.4)-(2.6)  is  written: 

dD    ~     2  ,- 

^  -  f C  -  H  gh  =  0  ,  (2.16) 

lr  +  f°  =  °  >  (2.17) 

dh 

ZT+HD  =  °  '  (2-18) 

2    2 
where  y   =  k  .   The  solution  to  this  set,  which  can  be  obtained  directly 

or  by  using  D  =  iku  and  £  =  ikv  in  (2 . 13) -(2 . 15) ,  is  given  by: 

fl                           y2gho 
D(k,t)  =  D   cos  vt  +  sin  vt  +  sin  Vt  ,       (2.19) 

~  v 

fD  2      2 

C(k,t) -j2-  sin  vt  +[1L-f-  +  -o  cos  vtJ  ^0 

V     v 
2 
-  -^-1  [1  -  cos  vt]  ghQ  ,         (2.20) 
v 
HD  f 

h(k,t)  = sin  Vt r  [1  -  cos  Vt]  C, 

V 

.2    2 


+l\  +  1LfH"  COS  V0  ho  '        (2-21) 


V      V 

where:  v  =  f  +  y  gH  .  (2.22) 


3.   Finite  Difference  and  Finite  Element  Solutions 

Schoenstadt  (1980)  carried  out  a  general  analysis  of  the  solutions  to 
(2.1)-(2.3)  which  allowed  for  spatially  centered  finite  differences  or 
finite  elements.   We  will  use  the  same  method  to  compare  certain  finite 
difference  and  finite  element  solutions  to  systems  (2.1)-(2.3)  and  (2.4)- 
(2.6).   The  various  finite  difference  and  finite  element  forms  correspond- 
ing to  (2.1)-(2.3)  or  (2.4)-(2.6)  are  given  in  the  Appendix.   Following 
Schoenstadt  (1980)  the  Fourier  transformed  versions  of  the  various 
numerical  schemes  for  the  equations  (2.1)-(2.3)  can  be  written  in  the 
following  form: 
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a(k)  ^r  =  fB(k)  v  -  iga(k)  h  ,  (3.1) 

at 

a(k)  ^7  =  -f3(k)  u  ,  (3.2) 

at 

a(k)  ^-  =   -iHa(k)  u  .  (3.3) 

at 


The  functions  a(k),  3(k)  and  o(k)  are  given  in  Table  I  for  the  various 
schemes  considered.   This  set  can  be  put  in  the  same  form  as  (2.8)-(2.10) 
by  dividing  by  a  and  by  setting: 

n  =  3/a   and   u  =  a/a  .  (3.4) 

In   this   case  the   frequency  is  given  by 

v2  =   (32f2  +  a2gH)/a2    .  (3.5) 

The   solutions   to   set   (3.D-C3.3)   are  given   by  (2.  12)-(2.  14)   with  the  use 
of   (3.4)    and    (3.5). 


Table   I.      Coefficients   in   primitive  equations   for   various 
numerical   schemes. 

Scheme  a 

differential  1 

A  1 

B  1 

FEM   A  (2+cos(kAx))/3        (2+cos(kAx) )/3 

FEM  B  (2+cos(kAx))/3        (2+cos(kAx ) )/3    (5sin(kAx/2)+sin(3kAx/2) )/4Ax 


6 

a 

1 

k 

1 

sin    (kAx)/Ax 

1 

sin    (kAx/2)/(Ax/2) 

(kAx))/3 

sin    (kAx)/Ax 

11 


The  numerical  schemes  for  the  vorticity-divergence  system  (2.4)-(2.6) 


lead  to  the  following  transformed  equations: 

a  —  -  af  £  -  a  gh  =  0  , 
at 


a  4r  +  fa°  =  °  » 
at 


(3.6) 
(3.7) 


dh 


a  —  +  HaD  -  0  , 
dt 


(3.8) 

where  a(k)  and  a(k)  are  given  in  Table  II.   This  set  can  be  put  in  the 
same  form  as  (2. 16 )— ( 2. 18)  by  dividing  by  a  and  setting: 


2    2. 

u  =  a  /a 


(3.9) 


The  frequency  equation  (2.22)  becomes 


v2  =  f2  +  (a2/a)gH  , 


(3.10) 


which  has  a  different  form  from  (3.5).   The  solutions  to  set  (3.6)- (3.8) 
are  given  by  (3.6)-(3.8)  with  the  use  of  (3.9)  and  (3.10). 


Table  II.   Coefficients  in  vorticity-divergence 

equations  for  various  numerical  schemes. 


Scheme 

Scheme 

differential 


Q' 


a 
1 
finite  difference      1 

FEM  (2+cos(kAx))/3     sin2(kAx/2)/( Ax/2) 2 


sin2(kAx/2)/(Ax/2)2 
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The  various  parameters  which  determine  the  solutions  (2. 13)— (2. 15)  and 
(2. 19 )- (2.21)  are  shown  in  Tables  I  and  II,  respectively.  Table  I  contains 
Schemes  A  and  B  for  the  primitive  equations  where  Scheme  A  is  unstaggered 
and  Scheme  B  has  the  velocity  points  midway  between  the  height  points  (see 
Schoenstadt,  1980).  The  table  also  includes  the  finite  element  forms 
which  are  obtained  when  piecewise  linear  basis  functions  are  used.   Note 
that  k  is  poorly  represented  by  cj  with  Scheme  A  near  k  =  tt/Ax,  and  that 
the  problem  remains  with  the  FEM  version  of  Scheme  A.   The  staggered  grid 
gives  a  much  better  approximation  since  spatial  derivatives  are  computed 
over  a  distance  of  Ax  compared  to  2Ax  with  the  unstaggered  grid. 

Table  II  contains  the  parameters  for  the  finite  difference  and  finite 
element  versions  of  the  vorticity-divergence  set  of  equations.   In  this 

case  vorticity,  divergence  and  height  are  carried  at  the  same  points. 

2  2 

Note  that  o  for  both  cases  is  the  same  as  the  value  of  a   for  Scheme  B 

from  Table  I.   It  can  be  seen  from  the  tables  that  the  staggered  primitive 

equation  (Scheme  B)  and  vorticity-divergence  formulations  have  the  same 

values  for  a  and  a  and  therefore  for  v,  so  that  these  should  give  the  same 

solution  except  for  truncation  error  in  the  initial  conditions. 

As  pointed  out  by  Schoenstadt  (1980),  the  solutions  (2. 12)-(2. 14 )  for 

the  various  schemes  differ  only  through  the  coefficients  lA*,  PA,  and 

2 
ny/v  ,  and  the  same  dependence  occurs  in  system  (2. 1 9 )— (2 . 21 )  with  n  =  1, 

except  that  the  coefficient  nu/v2  does  not  appear.   Figure  la  contains 
the  phase  velocity,   c  =  v/k  ,  as  a  function  of  kAx/fr   for  the  various 
schemes  in  Tables  I  and  II  as  computed  from  (3.5)  and  (3.10),  respec- 
tively.  The  differential  solution  approaches   f/k   for  small  k  and 
the  shallow-water  speed  (gH)1/2  for  large  k.   Scheme  A  gives  the  poorest 
phase  speed  and  the  finite  element  Scheme  A  is  also  very  poor  for  the 
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highest  wavenumbers.   The  finite  element  scheme  B  is  very  close  to  the 
differential  solution,  while  the  vorticity-divergence  FEM  scheme  is  a 
little  higher.  The  group  velocity,  G  =  dv/du,  is  given  in  Fig.  1b,  as  a 

function  of  kAx/ir.   The  differential  solution  is  zero  at  k  =  0  and  it 

1/2 
approaches  the  shallow-water  phase  speed  (gH)    for  large  k.   Scheme  A 

and  its  FEM  version  are  very  poor  for  the  short  waves  since  they 

propagate  energy  in  the  wrong  direction.   In  fact  the  FEM  scheme  gives  a 

group  velocity  which  is  more  than  double  the  correct  value  and  of  the 

wrong  sign,  at  certain  points.  The  FEM  scheme  B  gives  the  best  group 

velocity  while  the  FEM  vorticity-divergence  scheme  gives  values  that  are 

somewhat  higher. 

2 
The  coefficients  n/v,  u/v  and  nu/v  are  given  in  Fig.  2a,  2b  and  2c, 

respectively,  as  functions  of  kAx/iT.   Scheme  A  is  the  poorest  for  each 

coefficient,  but  the  FEM  version  of  scheme  A  is  just  as  bad  for  the  short 

waves.   The  best  scheme  is  the  FEM  version  of  scheme  B,  although  the  FEM 

vorticity-divergence  scheme  is  also  very  good.  The  coefficient  n/v,  which 

is  given  in  Fig.  2a,  is  especially  important  since  n/v  relates  the 

initial  height  to  the  final  (steady-state)  height  field  (see  (2.14)).   In 

particular,  the  figure  shows  that  if  v   =0,  the  final  h  for  k  =  tt/Ax  is 

more  than  25  times  too  large  for  scheme  A  and  the  FEM  version  of  A!   This 

is  one  reason  why  non-staggered  schemes  tend  to  generate  small  scale 

noise.  These  results  were  given  by  Schoenstadt  (1980)  with  the  exception 

of  the  vorticity-divergence  schemes. 

4.   Final  State  Example 

The  two  aspects  of  the  geostrophic  adjustment  process  that  must  be 

considered  in  assessing  a  particular  numerical  scheme  are:  1)  forecast 

time  required  to  reach  the  adjusted  state,  2)  the  accuracy  of  the  final 
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adjusted  state.   The  group  velocity  curves  in  Fig.  1a  provide  an 
indication  of  the  comparative  adjustment  times  for  the  various  schemes. 
The  final  adjusted  state,  which  is  more  important,  could  be  obtained  by 
Fourier  transforming  the  terms  that  are  independent  of  t  in  (2. 12)-(2. 14) 
or  (2. 1 9 )— (2 . 21 ) .   However,  in  this  paper  the  final  state  will  be 
determined  by  integrating  the  finite  difference  equations  in  t  until  the 
adjusted  state  is  reached.   This  approach  is  preferable  because  time 
differencing  effects  are  included  and  a  time  filter  can  also  be  used. 

The  various  sets  of  equations,  which  are  given  in  the  Appendix  are 
integrated  with  centered  time  differences.   The  time  filter  developed  by 
Robert  (1966)  (see  also  Asselin,  1972)  is  applied  to  the  past  time  value 
with  the  coefficient  y  =   .05.   The  new  time  values  for  the  FEM  schemes 
are  found  by  Gauss  elimination. 

The  initial  conditions  are  given  by: 

a   |x|  <  Ax/2 
h(x,o)  =  (  (4.1) 

o   |x|  >  Ax/2  , 

u(x,o)  =  v(x,o)  =  0  ,    or    c;(x,o)  =  D(x,o)  =  0  . 

These  initial  conditions  are  convenient  for  comparing  the  various  schemes 
since  no  truncation  error  is  introduced  when  the  initial  vorticity  and 
divergence  are  computed  from  these  initial  velocities.  The  analytic 
solution  for  the  final  adjusted  h  field  can  be  obtained  by  integrating  the 
following  expression  that  was  obtained  by  Schoenstadt  (1977): 


h  (x)  =h(x,o)  ."-A-   ;'   sgn(x-?)e-|x-?l/A[f  |^,0)-v(^,0)  ]d?  , 

Tl^r  I  OX 


/ 
_oo 

(4.2) 
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1/2 
where  h  (x)  is  the  final  adjusted  height  and  A  =  (gH)  '  /f  is  the  Rossby 

5 

radius  of  deformation.   The  initial  geostrophic  wind  which  is  required  in 
(4.2)  can  be  conveniently  written: 

§|^(x,0)  =  ^  [6(x+Ax/2)-6(x-Ax/2)]  ,     (4.3) 

f  dx  f 

where  <5(x)  is  the  delta  function. 

When  (4.1)  and  (4.3)  are  introduced  into  (4.2)  the  solution  becomes: 

e"x/A  sinh(Ax/2A)      Ax/2  <  x 

h  (x)  =  a  1  -  e_Ax/2A  cosh(xM)  -Ax/2  <  x  <  Ax/2  .    (4.4) 

s  —   — 

ex/A   sinh(Ax/2A)         x  <  -Ax/2 
Fig.  3  contains  h  (x)  for  the  case  Ax  =  A/2. 

The  numerical  integrations  with  the  various  schemes  are  performed  on  a 
grid  of  200  points  with  cyclic  boundary  conditions.   The  initial  distur- 
bance at  x  =  0  is  placed  in  the  center  of  the  computational  domain  so  that 
the  cyclic  boundary  conditions  will  not  affect  the  solution  near  x  =  0 
until  well  after  the  adjusted  state  is  reached.  Fig.  3  includes  the  numer- 
ical solutions  at  t  =  3  days  for  the  following  schemes:   A,  B  and  FEM  A. 
Scheme  A  shows  strong  oscillations  with  every  other  point  returning  to  0. 
The  FEM  scheme  A  has  smaller  oscillations  near  x  -  0,  but  they  become 
larger  than  the  oscillations  with  scheme  A  farther  out.   Scheme  B  gives 
very  smooth  behavior  and  it  is  close  to  the  analytic  solution.   The 
vorticity-divergence  system  gives  the  same  solution  as  scheme  B,  and  is 
very  close  to  the  analytic  solution  as  can  be  seen  in  Table  III  which 
compares  the  solutions  at  the  first  two  grid  points. 
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Table  III.   Numerical  solutions  h/a  at  t  =  72  hours  for  the  first  two  grid 
points  for  various  schemes  compared  with  analytic  solution. 


X 

0 

Ax 

Differential 

0.221 

0.153 

A 

0.459 

0.0 

B 

0.240 

0.148 

vorticity-divergence 

0.240 

0.  148 

FEM  A 

0.298 

0.084 

FEM  B 

0.227 

0.157 

FEM  vorticity-divergence 

0.213 

0.154 

The  results  given  in  Fig.  3  and  Table  III  are  consistent  with  the 

2  2 
curves  for  H /v  shown  in  Fig.  2a,  since  h  is  proportional  to  n  /v   (see 

s 

(2.14)  and  (2.21)).   In  particular  the  poor  behavior  for  the  unstaggered 
primitive  equation  schemes  (A  and  FEM  A)  in  Fig.  2a  is  consistent  with  the 
large  amplitude  short  waves  in  Fig.  3.   Also  the  large  oscillations 
farther  out  with  FEM  A  may  be  the  result  of  the  large  spurious  group 
velocity  that  is  shown  in  Fig.  1b  for  that  scheme.  All  the  staggered 
primitive  equation  and  vorticity-divergence  schemes  give  excellent 
predictions  of  the  final  adjusted  height  field.   It  should  be  pointed  out 
that  the  inclusion  of  light  time  smoothing  CY  =  .05)  is  necessary  to 
produce  the  spatially  smooth  solutions  for  these  cases.  Apparently  the 
vanishing  group  velocity  for  kAx/TT  =  1  (see  Fig.  1b)  does  not  allow  the 
smallest  scale  gravity  waves  to  propagate  out  from  the  initial  distur- 
bance.  Haltiner  and  McCollough  (1975)  demonstrated  the  usefulness  of  time 
filtering  in  a  baroclinic  primitive  equation  model. 
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5.   Conclusions 

The  objective  of  this  paper  is  to  determine  the  response  of  various 
finite  element  schemes  to  small  scale  initial  conditions  or  small  scale 
forcing.   It  is  especially  important  that  FEM  prediction  schemes  properly 
describe  small  scale  features,  because  FEM  models  usually  require  more 
computational  effort  per  degree  of  freedom  than  most  finite  difference 
models.   This  study  treated  the  geostrophic  adjustment  process  with  the 
linearized  primitive  equations  and  also  with  the  related  vorticity- 
divergence  set  of  equations.   The  development  followed  Schoenstadt  (1980) 
wherein  the  spatially  discretized  equations  were  Fourier  transformed  in  x, 
and  then  solved  with  arbitrary  initial  conditions.   These  solutions  were 
dependent  on  certain  coefficients  which  were  computed  for  the  various 
numerical  schemes  and  compared  with  the  differential  expressions.  Three 
FEM  schemes  were  examined  as  well  as  the  three  corresponding  finite 
difference  schemes.   It  was  found  that  the  unstaggered  (scheme  A)  primi- 
tive equation  model  gives  the  poorest  behavior  followed  by  the  correspond- 
ing FEM  formulation.   These  schemes  are  especially  bad  for  the  shortest 
resolvable  scales.  The  finite  difference  primitive  equation  model,  which 
staggers  height  points  between  velocity  points  (scheme  B)  has  much  better 
behavior  than  the  unstaggered  schemes.  The  vorticity-divergence  model 
where  C,  D  and  h  are  carried  at  the  same  points  has  the  same  coefficients 
as  scheme  B.  The  FEM  version  of  scheme  B,  which  has  staggered  nodal  points, 
was  found  to  have  the  best  behavior  and  the  FEM  vorticity-divergence  model 
was  also  found  to  be  very  good. 

The  six  schemes  were  also  compared  by  integrating  them  numerically  with 
centered  time  differences  from  an  initial  state  at  rest  with  a  height 
perturbation  at  a  single  point.   The  analytic  solution  for  this  initial 
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state  approached  a  smooth  height  field  after  the  inertial  gravity  waves 
radiated  away.  Scheme  A  and  the  FEM  form  of  scheme  A  gave  very  poor 
solutions  with  large  oscillations  from  point  to  point.  All  of  the  other 
schemes  produced  smooth  solutions  with  the  FEM  schemes  being  the  most 
accurate.   The  smoothness  of  these  solutions  was  improved  by  light  time 
smoothing.   Although  the  initial  state  used  in  this  comparison  is  somewhat 
extreme,  it  shows  clearly  the  superiority  of  the  staggered  primitive 
equation  and  vorticity-divergence  schemes  over  the  non-staggered  primitive 
equation  schemes. 

Winninghoff  (1968),  Arakawa  and  Lamb  (1977)  and  Schoenstadt  (1980) 
have  demonstrated  the  advantages  of  spatial  staggering  of  predictive 
variables  in  finite  difference  models.  Our  results  strongly  indicate  that 
FEM  models  should  either  use  staggered  nodal  points  in  the  primitive  equations 
or  unstaggered  nodal  points  in  the  vorticity-divergence  equations  (see 
also  Schoenstadt,  1980).   In  fact  Staniforth  and  Mitchell  (1977,  1978) 
have  developed  a  FEM  model  based  on  the  vorticity-divergence  form  of  the 
shallow-water  equations  that  produces  smooth  forecasts  with  only  time 
smoothing.   In  contrast,  Kelley  and  Williams  (1976)  obtained  very  noisy 
results  with  an  unstaggered  FEM  model  which  used  the  primitive  equations 
for  flow  in  a  channel.   If  non-staggered  finite  FEM  element  models  are 
used,  it  is  often  necessary  to  use  high  order  smoothing  to  damp  the  small 
scales  as  discussed  by  Cullen  (1976).  Thacker  (1978)  tested  a  finite 
element  formulation  of  the  linearized  shallow-water  equations  with 
staggered  nodal  points  and  he  obtained  smooth  solutions. 

Since  FEM  models  usually  require  more  computer  time  per  degree  of 
freedom,  it  is  very  important  for  the  numerical  scheme  used  to  be  accurate 
for  as  small  a  scale  as  possible.   In  this  paper  we  have  shown  that  the 
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usual  non-staggered  FEM  formulation  of  the  primitive  equations  gives  very 
poor  geostrophic  adjustment  for  small  scale  initial  conditions.   The  same 
conclusion  follows  for  small  scale  heating.  On  the  other  hand  either  the 
use  of  the  primitive  equations  with  staggered  nodal  points  or  the 
vorticity-divergence  equations  with  unstaggered  nodal  points  gives 
excellent  treatment  of  small  scale  features  in  the  geostrophic  adjustment 
process.  Clearly,  the  use  of  either  formulation  should  be  much  more 
efficient  than  the  unstaggered  primitive  equations,  even  when  the  latter 
have  smoothing  to  destroy  the  smallest  scale  features. 
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Appendix 

In  this  Appendix  the  spatially  discretized  prediction  equations  are 

given  for  each  scheme  with  x  =  mAx.   The  following  schemes  approximate  the 

primitive  equations  (2.1)-(2.3): 

Scheme  A 

9um   „     g(Vl  "  hm-l>    n 
W-  -   fVm  + 2A^ =  °  ' 

9v 

m  +  fu  =  0  , 


FEM  Scheme  A 


Scheme  B 


3t      m 


9h       u_j_i  ~  u  i 

m  +  H  (  mflOA  m"1)  =  0  . 


3t      v    2Ax 


M^-fMv  +    m+19.   m-1  =  0  , 
at  m         2Ax 

3v 

M  3-2  If  Mu  =  0  , 
dt         m 

3h  "_■!   _  u    .. 

m   m±  u  /  nrl"1    m-lN    _ 

M  "5T  +  H  < 2A^ >  =  °  • 


dt      m  Ax       "  u  » 


3v 

aT~  +  f  u  =  °  » 
ot      m 


3h 


_■  +  R(  tl/2  ^  m-l/2)  _  Q 
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FEM  Scheme  B 


M  3Um   f„     5rhnrfl/2   hm-l/2,  ,  3  ,hnr+l/2   hm-l/2. 

M  t: fMv  +  -r-( 7 )  +  t-7-( a )    °  > 

3t       m   8       Ax  16       Ax 

3v 

M  v^  +  f  Mu  =  0  , 
3t       m 

m  9hm  u  ur5  rUnH-l/2  "  Vl/2,    3  "m+3/2  "  "m-3/2,    . 
M  JT  +   H[8  ( Ax" >  +  8  ( Ax" >  =  °  ■ 


where  NP  =  (a    +  4a  +a  „)/6. 
m     m+1     m   m-1 

Scheme  B  is  staggered  in  such  a  way  that  the  height  points  are 
equi-distant  between  the  velocity  points.   The  FEM  equations  can  be 
derived  with  piecewise  linear  basis  functions  (see  for  example  Chapter  6 
in  Haltiner  and  Williams,  1980). 

The  vorticity-divergence  system  (2.2l)-(2.6)  is  approximated  with  the 

following  schemes. 

Vorti city-Divergence 

3D  h    .,    -   2h     +  h     . 

m        j._      .         /   m+l  m  m— 1*         _ 

JT  ~   f  Sn  +  8    ( 71 >    =   °    • 

Ax 

tj—  +  f D     =  0    , 

dt  m 

3h 

a  +  HD     =   0    . 


3t      m 
Finite  Element  Vorticity-Divergence 

3D  h,.    -   2h     +  h      , 

vt        m        ctjtr      j.  nn~l m  m-1   _   „ 

M  ^- fMC      +  g  5 0    , 

dt  m  A    / 

Ax 

Km 

M  tH11  +   f  MD      =   0    , 
dt  m 

3h 

M  7T-5  +  HMD     =   0    . 
dt  m 
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Figure  Captions 


Fig.  1.   The  phase  velocity  c  =  v/p,  and  the  group  velocity  G  =  dv/dy  as 
functions  of  kAx/iT  for  various  numerical  schemes.   The  curves 
are  labeled  as  follows:  1)  differential  solution,  2)  scheme  A, 

3)  scheme  B  and  vorticity  divergence  finite  difference  scheme, 

4)  FEM  scheme  A,  5)  FEM  scheme  B,  6)  FEM  vorticity-divergence 
scheme..  These  results  use  the  following  values: 

gH  =  10  m2s-2,  f  -  10_4s-1,  Ax  =  500  km. 

2 
Fig.  2.   The  coefficients  n/v,  u/v  and  ny/v  as  functions  cf  kAx/n,  with 

the  same  labeling  as  in  Fig.  1. 


Fig.  3.   The  numerical  solutions  for  schemes  A,  B  and  FEM  A  as  functions 
of  x/Ax  at  t  =  3  days.   The  steady-state  differential  solution, 
which  is  given  by  (4.4),  is  included  for  comparison. 
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